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Abstract 

This paper constructs perfectly matched layers (PML) for a system of 2D Coupled Nonlinear Schrodinger 
equations with mixed derivatives which arises in the modeling of gap solitons in nonlinear periodic structures 
with a non-separable linear part. The PML construction is performed in Laplace Fourier space via a modal 
analysis and can be viewed as a complex change of variables. The mixed derivatives cause the presence of 
waves with opposite phase and group velocities, which has previously been shown to cause instability of layer 
equations in certain types of hyperbolic problems. Nevertheless, here the PML is stable if the absorption 
function a lies below a specified threshold. The PML construction and analysis are carried out for the linear 
part of the system. Numerical tests are then performed in both the linear and nonlinear regimes checking 
convergence of the error with respect to the layer width and showing that the PML performs well even in 
many nonlinear simulations. 

Keywords: perfectly matched layers, coupled nonlinear Schrodinger equations, mixed derivatives, group velocity, 
stability 

1 Introduction 

Perfectly matched layers (PML) are a relatively simple and efficient tool for the truncation of spatial domains 
for wave type problems posed on unbounded (or large) domains. In numerical simulations such a truncation is 
often necessary as well as desired and PML guarantees that waves traveling through the boundary are absorbed 
and reflections that are only exponentially small with respect to the layer width occur. PML have been first 
proposed by Berenger [T] for Maxwell's equations and since then derived and analyzed for many other equations, 
like wave and Helmholtz equations [5] , linearized Euler equations [31 H] , general first order hyperbolic systems 
[S], Schrodinger equation [SJIT], etc. This paper proposes, analyzes and tests PML for a 2D Coupled Nonlinear 
Schrodinger system (CNLS) with mixed derivatives 

idtu, + {afdl + afdl + I3,d.,dy)u, + TU,{ui, . . . , txjv) = 0, j G {1, . . . , iV}, (LI) 

where A/} is a polynomial (typically cubic or cubic-quintic) nonlinearity and ctf \ '-^'f^^ l^i ^'^^ ^ '"'^^^ numbers. 
In [5], where it is called a system of Coupled Mode Equations, this system is shown to be an asymptotic model 
for gap solitons in the 2D periodic nonlinear Schrodinger equation with a finite contrast non-separable periodic 
potential. Previously the author together with T. Hagstrom have studied in [5] PML for ID and 2D Coupled 
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Mode Equations governing gap solitons in periodic structures with infinitesimal contrast. In that work the modal 
analysis in Laplace- Fourier space was used. The ID problem was hyperbolic and the general PML construction 
for first order hyperbolic systems [5] employing auxiliary variables was used. The 2D case was of a mixed type 
and required a combination of the method in [5] and a complex coordinate stretching. The problem ( 1 . 1 1 at 
hand is of generalized Schrodinger type and the modal analysis in Laplace-Fourier space reveals that a complex 
change of coordinates is sufficient for PML construction. In this paper (1.1) is studied under the condition 
a^^^'aj*'-* > /3| Vj e {1, . . . , N}, which is dictated by the asymptotic derivation in [8 and implies ellipticity of 
the spatial linear operator in ( 1 . 1 1 . 

As ( 1.1 1 is a model for nonlinear solitary waves, the following scenarios are particularly relevant for numerical 
investigations: evolution of a perturbed solitary wave; interaction of a solitary wave with a defect in the medium; 
or collision of several solitary waves. All these processes will typically lead to shedding of radiation that usually 
has small amplitude compared to the pulse(s), travels away from them, and needs to be treated at the boundary of 
the domain of interest. Examples of studies using PML in such situations are [inillTlIIl]. In addition, simulations 
where a pulse of magnitude comparable to the solution maximum leaves the domain are often desired. In such a 
case the polynomial nonlinearity cannot be in general neglected in the layers and leads to truly nonlinear layer 
dynamics. 

The derivation and analysis of PML in this paper is based purely on the linear part of the system 
nevertheless the presented numerical tests demonstrate satisfactory functionality even in prototypical examples 
corresponding to all of the above nonlinear scenarios. The analysis guarantees that in the linear regime the layer 
is absorbing and perfectly matched. Stability of the layer equations in time is shown to hold if the maximum of 
the absorption function a lies below a threshold, which diverges to infinity for /3 — s- 0. The layer equations are, 
therefore, conditionally stable, which is in spite of the presence of plane waves in the linear part of (1.1 1 with 
opposite group and phase velocities. Such a mismatch of group and phase velocities has been shown in certain 
hyperbolic systems to lead to instability of the layer equations [TSJ [Hj . PML for 3D linear Schrodinger equations 
with mixed derivatives have been previously used in |15j . Perfect matching and stability were, however, not 
analyzed there in the presence of mixed derivatives. 

The rest of the paper is organized as follows. In Section [2] relevance of the CNLS system is discussed and 
conditions on coefficients are provided which allow removal of the mixed derivatives via a change of variables. 
The PML is derived via the modal analysis in Laplace- Fourier space in Section |3] In Section |4] stability of the 
linear (F = 0) layer equations is analyzed. Finally, Section [5] presents a number of numerical tests in both the 
linear and nonlinear regimes. Exponentially fast convergence of the error within the physical domain with respect 
to the layer width is verified via linear tests and observed, though with larger error values, even in nonlinear 
tests. 



2 Relevance of the CNLS with Mixed Derivatives 



Systems of the type (1.1 1 have been shown in [8 to describe gap solitons in "nonscparable" Kerr nonlinear 
structures for values of their spectral parameter (frequency or propagation constant) lying in an asymptotic 
neighborhood of a spectral gap edge. The CNLS system is then called the Coupled Mode Equations (CMEs) 
and governs the dynamics of slowly varying envelopes of the gap soliton. The particular model for which CMEs 
were derived in [8 was the periodic nonlinear Schrodinger equation 

idt,t^ + AV^ - V{x', yyj + FIV-PV = 0, V{x' + d,,y') = y' + d^) = V{x\ y') V(a;', y') e (2.1) 

with some di^2 > and where the periodic structure V is fixed, i.e., does not depend on the asymptotic parameter, 
and non-separable. In physics literature this case is typically referred to as a large contrast periodic structure. 



Equation (2.1 ) describes propagation of light in 2D photonic crystals as well as evolution of matter waves in Bose 
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Einstein condensates. The nonseparability condition requires that V{x',y') ^ Vi(a;') + V^iy') for any functions 
Vi , V2 . For a gap soliton near a gap edge defined by N maxima or minima of the corresponding band structure 
u;(fc) the general CMEs read as with the coefficients OLf^^ oif^ and jij proportional to the second derivatives 
of the spectral bands at the extrema with respect to the components of the wavevector fc, see [51. In particular, 
fij is the mixed second derivative. As the extrema are either all minima or all maxima, the coefficients satisfy 
sign(aj'=^) = sign(aj^^) Vj G {1, . . . , A^} and 



afaf > P\ 



Vje{i,...,iV}, 



(2.2) 



which also guarantees ellipticity of the spatial differential operator in 

In contrast, structures with a separable linear part, like the periodic nonlinear Schrodinger equation (2.11 
with V{x,y) = Vi{x) + ¥2(1/), which was studied in [El [17], lead to f3j = Vj and the CMEs take the form of 
classical CNLS systems. 



CMEs of the type (1.1 1 can also be derived as an approximative model in the same asymptotic regime as 
above for gap solitons in the Maxwell problem Aip — V{x' ,y')d'^,ip — Tdf,{ip'^) = 0, ip{x,t) G M, with a finite 
contrast periodic structure V{x' , y') as these do not have a separable linear part either. 

In [8] an example of the potential V{x' , y') is presented, for which the band structure indeed leads to (3j 7^ 0. 



A simple prototypical example of the system (1.1 1, which is used in this paper for some of the numerical 
tests, is 



idtU2 + {a^2 dl + "2 <9y + P2d^^dy)u2 + T [|u2pM2 + (2|uipU2 + ulu2) + eq\u2tu2] =0, 



(2.3) 



where the quintic nonlinearity with < has been included in order to avoid issues with blowup of solutions 
of the cubic 2D nonlinear Schrodinger equation [TS] . 



Removal of the mixed derivative via a change of variables The mixed derivatives in (l.ll (and in (2.3l) 
can be in certain cases removed via the change of variables 



x\ j a cos 9 — 6sin6'\ jx 
A) I \asiiiO bcos9 I \y j 



(2.4) 



for some a,b S K, leading to 



cos^ + a^f &2 gj^2 e-p^"^ sin 26 ] dj 



a) a 



a'f^a^ sin^ 9 + af'b' cos' + (3j^ sin 261 ) ft 



ah 



+ ((af - a'f h^) sin 26* + /3ja6cos26l) didy. 



The mixed derivative is, clearly, removed if ~ ^^j-r = . . . = '^jtjt by the choice a — b i ^^-r ] and — ^■ 
Otherwise, the removal is successful if there are constants a, 6 € M such that 



1/2 



P2 



via the choice ^ — tan ^ 



/3j ah 



(2.5) 
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Note that in the case N = 2 the condition (2.5) reduces to 



4!l 



(2.6) 



which is solvable always unless ^ 



-Sign ( ^ - ^ 



(v) (y) / (^) 

and ^'^i ^^|— ^ or vice versa and unless sign ^^1— 

PI P2 ' ° \ PI P2 



As the mixed derivative in (1.1 1 cannot be removed in all cases, it is important to study PML for this system 
with Pj ^0. 



3 PML Derivation 

Since the derivation of perfectly matched layers and their analysis are performed only for the linear part of the 



CNLS system ( 1 . 1 1 , the analysis will be using merely the linear part of one scalar equation due to the fact that 



the system is diagonal in its linear part. The linear problem at hand, thus, reads 

idtu + {a^^'^dl + a^y^dl + I3d^dy)u = 0, (x, y) £ M^ t > 



(3.1) 



sign(Q!^^'') ~ sign(a*^^^). Note that the mixed derivative is not removed in (|3.1[) because the 



constructed PML will be used in the coupled system ( 1 . 1 1 , where the removal is not always possible as discussed 
in the previous section. 

Using the same approach as in [3] , the PML is first derived in the directions of coordinate axes x and y and 
then combining the x and y— layers, the corner layers are then proposed so that the resulting layer equations are 
applicable for all layers around the rectangular domain O as sketched in Fig. [l] 



Q 

physical domain 



L L +5 



Figure 1: Physical domain Q, ~ [Oji^:] x [0,Lj^] surrounded by layers. 

Without loss of generality the PML derivation is performed for the x— layer. For that end let us suppose the 
domain is unbounded in y, so that f2 = [0,La;] x IR- The Laplace transform in t with Re(s) > and Fourier 



transform in y (with the dual variable ky) of (3.1 1 yield 

(is + a^'^'^dl - a^y^l + i(3kyd^)u = 0. (3.2) 
The initial data in the Laplace transform vanish because of the assumption u{t = 0) = within the layers. The 
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modal solutions of (3.2 1 are 

u(a;;fc^,s) = e^^ A = A1.2 - ^ (-i/3fc^ ± ^-p^kl - 4a(-)(i.s - a(f)fc2) 
The ranges of A1.2 for Re(s) > are plotted in Fig. [2j 



(3.3) 





lm(l) 




-5' " t 











Figure 2: Ranges of Ai 2 for the modal solutions (3.3) 



The analysis using Laplace transform is helpful as it allows for immediate identification of modes with positive 
and negative group velocity Vg. Indeed, for Re(s) > propagating modes with > are contained within the 
modal set contained wholly in Re(A) < 0, i.e., for the problem at hand in the A2— set; and those with Vg <Q are 
contained within the modal set contained wholly in Re(A) > 0, i.e., in the Ai— set. This can be seen by performing 
a Fourier transform in x and Taylor expanding the dispersion uj{kx) about /cq G M for modes g'C^^"^^"*), where 
we set ikx — A and —ilu — s 



For (3.1 1 the above relation between sign(Re(A)) and sign(i;g) can be easily checked by studying the dispersion 
relation explicitly. For the propagating modes e'^'^="^~'^*-' with kx, w G M the relation reads uj = a^^^'k'^ 



^y'^ki 



Pkxky so that Vg{kx) — 2a'-^^kx + (iky. Clearly, > if and only if k^ > 



kx € (0, 2a5) ) i*^^ G ( 2cfw ' when > 0) the group velocity and the phase velocity Vp(kx) = uj{kx)/kx 

have opposite signsl This interval corresponds to the segment of the imaginary axis in Fig. [2] between the origin 



I3k 



-I3ky 



- a 



Let us stress that in the interval 



and the point 



since A = \kx- For many hyperbolic systems, like linearized Euler equations or equations 



for elastic waves, equality of sign(wp) and sign(i'g) is shown to be a necessary condition for stability of the PML 
equations in time [T31[T3]. Section |4] shows that in the Schrodinger type problem at hand this is not the case 
and stability can be easily achieved by a choice of the PML parameters. 



Let us now return to the transformed problem (3.2 1. For absorption in the a;— layers x < and x > Lx the 
solution ii needs to be modified (in the layers) to yield some u^^'^^ which in the layers satisfies 



dxU^ 



with 



Re(A) < for modes traveling right in x 
Re(A) > for modes traveling left in x. 



(3.4) 



As the A in (|3.3 1 satisfy the non-strict version of these inequalities, it would suffice to apply a simple rotation of A 
i^ by an angle p G (0, 7r/2), i.e., A = e'" (; 



about Ao 



The resulting modes ii^^^ — 



are, however, not perfectly matched with u at the interfaces x 
matching one can, instead, set 



and x = Lx- In order to achieve perfect 



-PML 



(3.5) 
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where Xq = and Xq = for the layers x < and x > L^c respectively. Choosing ^^cr^(a;o) — for all fc < n — 1 



guarantees C" matching of u and ii^'^^ a.t x = xq. With (3.5) one obtains A = A + ^A + i ^^^") j e^i^axix) in 
(3.4|. Under the condition rJx{x) > for x < and x > this A can be seen using (3.3l or Fig. b^to satisfy 



the inequahties in (p^If . Note that u^^^ 

0. 

,PML Yet US express dxU in terms of ii 



can be viewed as the solution over the whole domain [~5x, 



5 J if 



for x € [0,Lx] one sets crx{x) 
To derive equations for 



1 



where 



A + i 



pip 



/X, crxiOdS.- Thus, defining 9^ 



common factor e^'^, the x— layer equation reads 



{dx 



and eliminating the 



(3.6) 



/3 



2a(!') '''a^a: 



Treatment of the y— layers is completely analogous and defines the operator 9™"" .— i^^'Pa 

with ay ^ for y G [0, Ly],ay{y) > for ?/ < and y > Ly and with the perfect matching condition ^ay{yo) — 
for all fc < n — 1 at yo = and yo = Ly with a chosen n e N. The most natural approach to the corner layers 
[-5x,0) X [Sy,0), [-Sx,0) X {Ly,Ly + Sy], {L ^ + 5 ^ , L x] X (-5y,0) and {Lx + 5x,Lx] x (iy,L.y + 5j^] is to combine 
the two layer equations into 



= 0. 



(3.7) 



Note that as 



(a 



iPML^PML ^PML^PMLN ^^P^L _ C^^^ [3 

y y / 2 



the operators d™^ and 9,™''' do not commute unless jS — or (Jx,y are constant. Although in the performed 
numerical examples this does not seem to affect the L? error of the solution, in the numerical examples presented 
in Section |5] the operator 9™^ 9™^ was replaced by the commuting alternative \ (9™^9™^ + 9™^9™^). 

In the nonlinear case F 7^ the above analysis does not, strictly speaking, apply. Nevertheless, if the solution 
remains small within the layers, the polynomial nonlinearity can be viewed as negligible and it makes sense to 
simply use 



a 



^^^(97^)2 + /3,9 



PML^PML 



'■N 



^) = 0, J e {!,..., N} (3.8) 



as the PMF system corresponding to ( |F1| in such a nonlinear scenario. This formulation is used in the nonlinear 
numerical simulations in Section |5.2[ In fact, in one of the nonlinear numerical tests a large pulse enters the 
layer and the solution still qualitatively correct. 

To the author's knowledge no truly perfectly matched layers for nonlinear systems exist in the literature. 
Appending the linear layer equations with the corresponding nonlinear terms is a common approach. For the 
nonlinear Schrodinger equation (NLS) this is done, for instance, in [20]. In [2T] the same type of PML is 
constructed non-rigorously by viewing the nonlinearity as a (solution dependent) potential. It is then argued 
that the success of this approach is due to the time-transverse invariant property of NLS. Radiation boundary 
conditions, on the other hand, have been successfully derived for some truly nonlinear systems including NLS 
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4 Stability of the Layer Equations 



The analysis of Section [3] does not guarantee that the layer equations are stable in time. Stability is determined 
below only for constant Ux and ay, in which case rewriting the layer equations in Fourier space results in a 
diagonal ODE system. The following analysis determines which parameter values (in particular (Tx^y) lead to 
boundedness of all Fourier modes in time, i.e., stability, and which lead to growth of at least one mode, i.e., 



instability. Once again, as the linear system is uncoupled, it is sufficient to study the scalar problem (3.7 1 



4.1 Corner layer equations 

Let us firstly perform the change of variables x — |a*^^^ j^-'^/^a:, y — ja^^^ j^^/^y, which replaces a^^^ 1,0;'^^-' 1 
and P la'-^-'l^^/^ja'^^'l^^/^/? =: f3. Due to (2.2| one gets \(3\ < 1. Dropping the tildes over the spatial variables, 
one obtains for the corner layer equations 



(4.1) 



Applying the Fourier transform in x and y, with dual variables and ky, to (4.1 1 yields 



^ e'Pffx e'Poy 4 + e^'^pJi'^axOy 



v{kx, ky)if 



with /ia; = 1 + e^Pax and fiy — \ + e^Pffy. The stability requirement is Yle{v{kx, ky)) < y{kx, ky) E R-^. For the 
sake of simplicity let us set ax ^ ^y cr and p = 7r/4. Under these simplifications 



Re(i^) 



4(cr2 + 72cr+l)2 



{kl + kl) (2V2f3^a^ + aip^ - 4) - V2{P^ + 4) 
+kxkyP (V2cr^(/32 + 4) + a0^ - 4) - 8%/2) 



Clearly, for /3 ^ taking a > small enough yields Ke{iy) < while large ct > result in Re(i^) > 0. Next, 
Re(j') = if and only if 



— "-!/,32/,„2 



2/32(2CT2-l) + v^cr(/32_4)_g 



(75 - '^'S + ± (^'-' + ^-(/^V^ - 4) + (f - 2) - 4) 



1/2- 



Clearly, nonexistence of real solutions {kx, ky) is equivalent to D :— P'^a^ + \f2a(JP'a'^ ^ 4) + — 2j — 4 < 0. 
The roots of D are 



or,m = ^ ( 2 - /3 ± + 12/3 + 4 ) , a3,4(/9) 



72 

4^ 



2 + /3± - 12/3 + 4 



(4.2) 



A straightforward analysis of cti^2,3,4 reveals that I? < holds for < /3 < 1 when < o" < cri(/3) and for 
— 1 < /? < when < cr < a^dS). Because cr3(/3) = ai{~~l3), this reduces to 



< o- < (Ti(/9) for 



< 1. 



(4.3) 
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Figure § plots the function cri(/3). 




Figure 3: The stabihty threshold function (Ji{f3) in (4.2 1 



a 



(y) 



As a conclusion, under the condition (4.3 1 the layer equations (3.7) with constant — (Jy — a, p = 7r/4, a 



1 and \P\ < 1 are thus stable. Note that (Ti(/3) ^ oo as /3 ^ so that the layer equations are 



unconditionally stable for = 0, i.e., for the classical 2D Schrodinger equation. 



In order to determine the stability condition for the linear system (3.8) with F = each equation can be 



first scaled so that the coefficients of the non-mixed derivatives become 1 and the mixed derivative coefficients 



become :— 



A^)\-i/2iAv) 



ctj \ ^^'^Pj- The stability condition for {ax)j = {<^y)j ^ cr = const, is then 



< (7 < CTi ( max \Bi 
\je{i,...,N}' 



(4.4) 



Because in practice a^^y are usually taken non-constant, the above condition on a translates to a condition on 
max((Ta;) and max(crj,). This can be justified by approximating ax,y by piecewise constant functions and applying 



the stability condition (4.4) on each piece 



Clearly, the presence of waves with group and phase velocities of opposite sign (see Section |3| does not lead to 
unconditional instability of the layer equations (3.7), which is in contrast with the studies of several hyperbolic 
systems plfTl]. 



4.2 Side layer equations 



Taking the Fourier transform of the side layer equation (3.6) under the assumption of 
defining v analogously to Section |4.1| gives 



a =const. and 



Re(i/) 



cr(\/2+a-) 



{2kx 



SO that Re{v) < \/{kx, ky) € M^. 

Side layers are, therefore, unconditionally stable even for /3 ^ and it is possible to use an absorption 



function ax, whose maximum in the corners satisfies (4.4) and takes a larger value in the side layers y € [0, Ly\ 
(and analogously for ay) so that the absorption in the side layers is strengthened. This is, however, not done in 
the simulations in Section |5] and ax and ay are kept y and x independent respectively. 



5 Numerical Tests 



This section presents results of several numerical simulations of the system (3.8) in both the linear (F — 0) 
and the nonlinear (F ^ 0) case. The primary objective is to demonstrate convergence of the solution error 
with respect to the layer width. In the linear regime layers of infinite width {5x — 5y — oo) do not generate 
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any error and the restrictions of the solution of (1.1 1 and the solution of (3.8) onto the physical domain — 
[0, L^] X [0, Ly] are identical. Finite layers produce reflections from their far end but the resulting error in O decays 
exponentially with the layer width [531 US due to the exponential decay of the solution within the layers. 
This exponential error convergence is numerically demonstrated to hold also here. Even in the nonlinear tests, 
where exponential convergence cannot be proved for the proposed layer equations, the resulting convergence is 
apparently exponential although the relative error in the example with large data in the layers becomes large. 



In all the presented numerical examples the rectangular domain 



'^xi Lx+Sx] X [—Sy, Ly+Sy] wlth 6 



Lx, Ly > 



was used with = [0, L^] x [0, Ly] being the physical domain and the rest being the PML layers. The spatial 
discretization of the PDEs was done via the centered 4th order finite difference formulas 



30uij + 16ui+ij - Ui+2j)/(12dx^), 



dxu{xi,yj) « {ui-2,j - 8Mi_ij + 8ui+ij - Uj+2j)/(12rfa;), 



(5.1) 



where dx = — Xi and u^j = u{xi,yj)\ and analogously for the y— derivatives. The zero Dirichlet boundary 
condition was imposed at the outer layer boundary. The time-evolution was approximated via 4th order additive 
Runge-Kutta scheme of the ESDIRK type [27], in which the linear (stiff) terms are treated implicitly and the 
nonlinear terms explicitly. 

Regarding the PML parameters, p was taken p = tt/A and the absorption functions cTx and ay were chosen 
of the form 



(Tx{x) 



^[1 + tanh(a^(4)(a; - - + tim\i{Qax{5x){x - 

'^[1 - tanh(a,(4)(x + |^))][1 - tanh(6a,(5,)(a; + %))] 



'-))] for X e {Lx,Lx + 
for X e [—Sx,0) 



(5.2) 



with ax{6x:) = 12/Sx and analogously for (Ty{y)- Clearly, maxux = ^[l + tanh(aa;((52:)^)]-[l+tanh(a2;((52:)^^))] 
is well approximated by hx even for moderate values of 5x- The plot of Gx for Lxj&x = 5 is in Fig. |4] Note 
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Figure 4: The auxiliary PML function (Tj.(a;) in (5.2) 



that the function Ox can certainly be chosen differently than (5.2) and no claim on optimality is made here. The 



function in the second pair of square brackets on each line of (5.2) is used merely to make Ox converge to at 
X = and a; = in a smoother manner. One could, of course, drop this function and simply make ax{5-x) larger 
but that would result in a large slope of Gx within the layer, which leads to reflections in the numerical solution. 
An optimization study on the PML parameters, primarily Ox and CTj,, can be performed [28 to increase efficiency 
of the layers. 



5.1 Simulations of the Linear Case F = 



Clearly, in the linear case F = the system (3.8) decouples and the change of variables (2.4) (possibly distinct 
for each f) can be applied to remove the cross-derivatives. Nevertheless, because in the to-be-studied nonlinear 
case this removal is not always possible (see Section |2]), one of the two linear numerical tests provided below is 



9 



with the mixed derivatives present. Initial data locahzed at the center of f2 were used in the tests. The error was 
computed at < = 1 using the reference solution u^^^ determined by numerically solving the initial value problem 
( 3.8 ) with iV = 1 , r = via Fourier transform on a much larger domain (than fl) on which the dispersed solution 



at i = 1 is well localized. 



5.1.1 Linear Scalar Case with /3 = 

In the linear scalar case F = 0, = 1 with /3 = the problem reduces to the 2D linear Schrodinger equation 
and the layer equations (3.7 1 are those used extensively in the literature, see e.g. [TJUHIIII]- The numerical 



test for this case is presented here for completeness and comparison with the case /3 ^ as well as with other 
publications. 

The remaining coefficients are chosen a^^^ — 3/4 and a^^^ =5/4 and the numerical parameters are = Ly = 
6, dx — dy — « 0.017 and dt = 0.01. The PML parameters are p — 7r/4, hr^ — 30 and the computations 

were performed for 6 different layer widths dx — 5y € {0.08Lx,0.12Lx,0-lQLx,0.2Lx,0-25Lx,0.3Lx} with the 
initial data u{x, y, 0) = e 



-{x-L^/2y-{y-Ly/2y 



. Fig. [i] shows convergence of the norm of the error over the 
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Figure 5: Error convergence with respect to the layer width for the test in Section 5.1.1 stars: relative error 
e, = Wu'"'^^ - M'^^1U2(o)/||u'^^1|L2(n) at t = 1, dashed line: c lO'i-^^ 

physical domain at i = 1 in dependence on Sx — Sy featuring an exponential convergence e~^^'' with p « 1.82. 
5.1.2 Linear Scalar Case with (3 = 0.5 



The same parameters and initial data as in Section 5.1.1 are chosen here except for the following: a^^^ = a^^^ = 



l,/3 = 0.5 and the maximum value of cr^.y, which is set to hx — hy ~ 3.3, i.e., close to the stability threshold 
cri(0.5) « 3.325 in Q. 

As one can see in Fig. [e] the error convergence is, once again, exponential like e with p « 1.07. Compared 
to the /3 = case in Section [5. 1.1| the convergence is slower and the error values are slightly larger which is to be 
expected due to the weaker applied absorption. Fig. [7] shows the initial profile and the modulus of the solution 
at t = 1 for 8x = dy = 0-2Lx = 1.2. A simulation with the same coefficients and similar discretization and PML 
parameters to those in Fig. Q but with hx = hy = 20, which exceeds the stability threshold (4.3 1, is shown in 
Fig. [8j It, indeed, results in an instability within the layers, clearly seen in the plot at t = 0.6. 

5.2 Simulations of the Nonlinear Case T ^ 



As advertised in Section |2] the CNLS system (2.3 1 was used for numerical tests in the nonlinear case and 



Eq was taken negative in order to prevent possible finite time blowup of the solution [18'. Three tests are 



performed below. In the first two tests (Sections 5.2.1 and 5.2.2 1 the following choice of coefficients and numerical 
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Figure 6: Error convergence with respect to the layer width for the test in Section 5.1.2 stars: relative error 
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(b) t = 1 
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Figure 7: The solution modulus for the test in Section 5.1.2 with 5x = = O-'^.L^. Black lines denote the 
interface between the physical domain and PML layers, (a) initial data; (b) solution modulus at i = 1; (c) 
solution modulus at t = 1, y = 4. 



(a) t = 0.4 



(b) t = 0.6 





Figure 8: The solution modulus for the unstable test in Section 5.1.2 
at i = 0.4; (b) solution modulus at t = 0.6. 



with = hy = 20. (a) solution modulus 



parameters was made: F = 0.5, Sq = — 0.2,La; = Ly ~ 14, dx = dy ^ 14/250 « 0.056, = 0.01, p = 7r/4 
and 5x = 5y & {O.OSLa;, 0.12La;, 0.16La:, 0.2ix, 0.25^^:, O.Sia,}. Section 5.2.1 presents a case where the choice of 
oi{\, af\ and /3i^2 allows a change of variables that removes the cross-derivatives and section 



5.2.2 



a case where 

this is impossible. In both cases the initial data are the sum of a stationary solitary wave and four perturbing 



11 



Gaussians so that the dynamics result in a large amount of radiation shed toward the boundary with a solitary 
waves remaining at the domain center. The solution is evolved up to i = 5. In detail, the initial data are 



t{x, y, 0) ^<l>s{x- LJ2, y - Ly/2) + 0. 



,y^g-2((x-p,)^ + (y-gfc)2) A 



(5.3) 



with pi 



2,3,4 



2 



2 



4 



3Lx 
4 



and gi, 2,3,4 = 



4 



2 



respectively, and where (f)g(x,y) is the positive 



spatial profile of the stationary solitary wave (ground state) e~^*'(t)s{x,y) of (2.3 1. (j)s{x,y) was computed via 
Newton's iteration on the corresponding stationary system, i.e., on (2.3 1 with idt replaced by 1 and with zero 
Dirichlet boundary conditions at 90. In detail (2.3) is first solved with (3i ^ (32 — for the radially symmetric 



Townes soliton with ui = U2 via the shooting method. Next, the solution is numerically continued via homotopy 
in (3i and /32 solving for ui and U2 via Newton's iteration at each /3-step and using the previous solution as an 
initial guess. 



The third example (Section 5.2.31 tests the designed PML for the scenario of a pulse entering the layers in 
the nonlinear regime; for the choice of parameters see the corresponding Section. 

For all three tests exponential convergence of the error is observed despite the fact that the problem is 
nonlinear and there is no guarantee for such a convergence. In the first two tests the relative error is satisfactory 
while in the third case which is truly nonlinear even in the layers, the relative error is rather large but the solution 
is still qualitatively correct. Note that for all three tests below the figures with the solution profiles show only 
the first component ui as U2 behaves in qualitatively the same way. 



5.2.1 Nonlinear System with (3j — 



(a:) _ Ay) _ 3 Ax) _ Ay) _ 



Jy) 



Jy) 



and (3i — (32 = 0, which can be viewed as obtained 



= 1 and (3i — (32 — h via the transformation (2.4) with a = b = 1 



The coefficients here are al 

(x) (x) 

from the system with a\ — a2 
and 9 = J. The magnitude of the absorption functions ax,y is hx = hy = 8. 

Fig. [9] presents the error convergence at i = 5 with respect to the layer width, where the solution computed 
with the widest layer {S^ = 5y = 0.3^2;) was used as the reference solution u^^^ . The convergence is exponential, 
like 6"^"^^ with p w 0.47. Fig. 10 shows for 5^ ^ 5y = 0.2Lx = 2.8 the initial data and the modulus of the first 
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Figure 9: Error convergence with respect to the layer width for the test in Section 5.2.1 stars: relative error 

/||,,REF|| 
L2(0)/ll" llL2(f2) 



= llw™!- - M'^^''||r2roJlk'"=''L2m^ at t = 5, dashed line: c 10-°-"^^ 



component ui at times t = 2, when a large amount of radiation is traveling into the layers and ai t — 5, when 
most radiation has been absorbed. 
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(a) t = 



(b) t = 2 



(c) t = 5 





Figure 10: The solution modulus for the test in Section 5.2.1 with 5^ — Sy ~ Q.2Lj.. Black lines denote the 
interface between the physical domain and PML layers, (a) initial data; (b) solution modulus at i = 2; (c) 
solution modulus at t = 5; (d) and (e) solution modulus along y = Lj,/2 at i = 2 and t = b respectively. 



5.2.2 Nonlinear System with (3j ^ 



This example presents the case ai 



(x) 



1, a. 



{x) 



4' ^1 



a. 



iv) 



1 and Pi = 0.2,(32 = 0.15, for which 
equation (|2.6|) cannot be solved for nonzero a, b and thus the mixed derivatives cannot be removed. The i nitial 

differs from that used in Section 



condition (5.3 1 is, once again, used, where the solitary wave profile 



5.2.1 



due to the different PDE coefficients. The stability threshold for the selected coeflicients as given by (4.4) is 
cri(max(0.2, 0.15^)) — cri(0.2) « 7.67. In the simulation = hy = 7.6 was used. 



Using, once again, the solution computed with the widest layer (5^ 
the error convergence is plotted in Fig. 
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O.SLx) as the reference solution, 
"'jP Ki 0.35 can be observed. The 



where exponential convergence 
convergence rate is slightly smaller than that in Section [5. 2. 1| due to the weaker applied absorption. Fig. [T2]then 
shows the solution modulus at selected instances of time. Fig. 12 shows for Sx — Sy — 0.2Lx — 2.8 the initial 
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Figure 11: Error convergence with respect to the layer width for the test in Section 5.2.2 stars: relative error 
= IIt.--- - u'^^^\\mn)/\\u''^^\\LHn) at t = 5, dashed line: c lO-^-^s . 
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data and the modulus of the first component ui at times t = 2, when a large amount of radiation is traveling 
into the layers and t = 5 when most radiation has been absorbed. 



(a) t = 



(b) t = 2 



(c) t = 5 
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Figure 12: The solution modulus for the test in Section 5.2.2 with 5x — Sy — 0.2Lx- Black lines denote the 
interface between the physical domain and PML layers, (a) initial data; (b) solution modulus dX t — 2\ (c) 
solution modulus at t = 5; (d) and (e) solution modulus along y = Lj,/2 at t = 2 and t = b respectively. 



Finally, to check the stability result in long time dynamics, Fig. [13] shows the solution profile sX t — 200 and 
demonstrates that no growth occurs. 

t = 200 




Figure 13: The solution from Fig. 12 at i = 200 



5.2.3 Nonlinear System with (3j ^ and a Pulse Propagating into the Layer 



In order to test the performance of the proposed layer equations in a truly nonlinear regime. Figs. 14 and 15 
present the case of a pulse propagating into the layer (entering it at a corner of 51). The PDE coefficients are 
taken the same as in Sec. |5.2.2| and the initial data used were the stationary solitary wave used in Sec. |5.2.2| 
centered at (a;,?/) — {Lx/2, Ly/2) and multiplied by a plane wave in order to induce motion of the pulse along 



14 



the y — X line 



uix, y, 0) = (j)s{x ~ y - Ly/2)e 



&i{(x-L.^/2) + {y-Ly/2)) 



(5.4) 



This initial condition does not correspond to an exact moving solitary wave solution of ( |2.3[ ) (such solutions have 
not been found for this system with [3i ^ P2) but the resulting solution propagates in a close to solitary manner. 

The physical domain f2 = [0,La;] x [0,Ly] was set by — Ly — 10, the spatial discretization by dx = dy = 
Lx/lSO « 0.056 and the simulation was carried out for six different layer widths 6x = Sy G {0.1Lx,0-15Lx, 0.2Lx, 



0.25Lx,0.3Lj:,0.35Lx}- The rest of the PML parameters was as in Section 5.2.2 Using the solution with 
6x = 0.35Lj; as the reference solution, the convergence of the relative error over is plotted in Fig. 14 at 
both t = 0.5 when the pulse has just entered the layer and at t = 3 when the pulse has propagated far from 
the domain and the exact solution in f2 consists only of an exponentially small tail plus small radiation due 
to the fact that the pulse is not an exact traveling wave. The relative error is much larger than in the previous 
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Figure 14: Error convergence with respect to the layer width for the test in Section |5.2.3| on the left error at 
t — 0.5 and on the right at t = 3. stars: relative -L^ error = — u^^^ \\ {n) / W (n) 1 dashed line: 

c 10^"-^^ ''^ and c 10""-^^ on the left and right respectively. 

examples, which were linear or effectively linear in the layers, but on the selected 6x range the convergence seems 



to be again exponential at both t = 0.5 and t — 3. More importantly, as one can see in Fig. 15 the qualitative 
behavior of the solution is correctly captured and the pulse, which has amplitude about 0.99 before entering the 
layers, leaves fl with only very small reflected waves (amplitude ^ 10^^) remaining. Fig. 15 shows the solution 



modulus at several instances of time over the whole spatial domain as well as along the line y — x, along which 
the pulse propagates. The first shown instance is at i = 0.5 because for t < 0.5 the pulse is simply traveling from 
its initial location at {Lx/2, Ly/2) toward the corner layer. 

The results of this test suggest that the layer equations can be applied even in many truly nonlinear cases, 
where the solution amplitude within the layers may become large, mainly if only qualitative behavior of the 
pulses is required. No guarantee can, however, be given that its performance will be satisfactory in all such 
cases. Examples of relevant nonlinear problems are interaction of several pulses or interaction of pulses with 
localized defects, where one or more pulses leaves f2 within the simulation time. 



6 Discussion 

The presented analysis of PML for the 2D Schrodinger equation with cross derivatives shows that the presence 
of the cross derivatives leads to the existence of linear (Fourier) modes with opposite group and phase velocities. 
Unlike in some hyperbolic systems [12 the resulting layer equations are only conditionally unstable and 
a choice of the damping functions (Jx,ay below a calculatable threshold leads to stability of the linear PML 
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Figure 15: The solution modulus for the test in Section 5.2.3 with 5x = Sy = 0.2Lj;. Black lines denote the 
interface between the physical domain and PML layers, (a) - (d) solution modulus at t = 0.5, 1, 1.5 and t — 3 
resp.; (e) - (h) solution modulus along y = x corresponding to (a)-(d) respectively. 



equations. Note that the damping of the PML is ensured based on an analysis of group velocity rather than 
phase velocity of the linear modes. 

In the nonlinear case the linear PML equations are simply appended with the (polynomial) nonlinear terms. 
The layer performance is then affected only slightly if the solution remains small in the layers as seen in the 
provided numerical tests. An analysis of the layer performance in the nonlinear case would be valuable. Of 
tremendous interest would a perfectly matched layer for truly nonlinear waves, i.e. without the smallness 
assumption. Such analysis does not appear in the literature. In the area of radiation boundary conditions, on 
the other hand, limited results for nonlinear equations exist, see [22l [23] . 

The paper studies PML in the 2D case. Nevertheless, in 3D the derivation is completely analogous. The 



linear scalar equation corresponding to (3.1 1 is in 3D 



(6.1) 



The modal solutions analogous to (3.3 1 are u{x; ky, fc^, s) — e with 



A = Ai,2 = [-liPiky + f32h) ± ^-{Piky + f32KY - 4a(-)(is - aiv^kl - a(-)k^, - P^kyk,)^ 

and, thus, 5™^ generalizes to 9™^ := [d^ - ^{f3idy + ^2dS)- The operators 9™^ and 9™^ are 

defined analogously and the layer equations are similarly to (3.7 1 and (3.8 1 obtained by replacing dx,dy and dz 
by 9™^ a^"^ and 9™^ respectively. 

The algebra in the stability analysis becomes in 3D, however, much more complicated and will be left for 
future work. Note that PML for 3D Schrodinger equations with mixed derivatives have been previously used in 
|15j . Perfect matching and stability were, however, not analyzed there in the presence of mixed derivatives. 
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